#include "Tools.h"

//constructor
DataTransfer::DataTransfer(pMesh mesh1, pMesh mesh2){

	int dim = mesh1->getDim();
	int columns= dim + 1;
	int rows = M_numVertices(mesh1) + 1;
	int size = M_numVertices(mesh2) + 1;

	//create array to allocate the node values
	pNodeValue = new double [rows];

	//create array to allocate the gradients
	pGrad = new double *[rows];
	for(int i=0;i<rows;i++) pGrad[i]= new double [5];

	//create matrix to allocate the geometric coefficients
	pGeomCoeff = new double *[size];
	for(int i=0;i<size;i++) pGeomCoeff[i]= new double [columns];

	//create array to allocate the linear and th quadratic interpolation result
	pInterpolatedVals= new double [size];

	
	  //create octree on the first mesh
  theOctree = OctreeCreate2<iterall>(mesh1->beginall(dim),mesh1->endall(dim),dim);
};

//destructor
DataTransfer::~DataTransfer(){};


void DataTransfer::fine_to_coarse(pMesh fine, pMesh coarse){

	//variables
	int dim = fine->getDim();

	switch (getRestrictionOperator()){

	case PURE_INJECTION: {
		cout << "PURE_INJECTION\n";
		calculate_GeomCoeff(coarse,dim);
		calculate_LinearInterpolation(coarse,fine,dim);
		break;
	}
	case HALF_WEIGHTING: {
		calculate_GeomCoeff(coarse,dim);
		calculate_LinearInterpolation(coarse,fine,dim);
		addNeighborsWeighting(coarse);
		break;
	}
	case FULL_WEIGHTING: {
		calculate_GeomCoeff(coarse,dim);
		calculate_LinearInterpolation(coarse,fine,dim);
		addNeighborsWeighting(coarse);
		break;
	}
	default:
		cout << "Restriction Operator unknown.\n";
		throw 1;
	}
}

void DataTransfer::coarse_to_fine(pMesh coarse, pMesh fine){

	//variables
	int dim = coarse->getDim();

	switch (getInterpolationOperator()){

	case LINEAR: {
		calculate_GeomCoeff(fine,dim);
		calculate_LinearInterpolation(fine,coarse,dim);
		break;
	}
	case QUADRATIC: {
		calculate_GeomCoeff(fine,dim);
		calculate_LinearInterpolation(fine,coarse,dim);
		calculate_Gradients(coarse);
		calculate_QuadraticInterpolation(fine,dim);
		break;
	}
	case ADAPTATIVE: {
		break;
	}
	case CONSERVATIVE: {
		break;
	}
	default:
		cout << "Interpolation Operator unknown.\n";
		throw 1;
	}
}


void DataTransfer::calculate_GeomCoeff(pMesh mesh, int dim) {

	if(dim==2){

		//variables
		pEntity vertice_1,vertice_2,vertice_3; //element vertices
		double ptn1[3],ptn2[3],ptn3[3]; //coordinates from the element nodes

		//Loop over vertices
		VIter vit = M_vertexIter(mesh);
		while ( pEntity vertex = VIter_next(vit) ){

			// get coordinate of a point from the mesh to be interpolated
			double xyz[3];
			V_coord(vertex,xyz);
			
			//cout<<xyz[0]<<" | "<<xyz[1]<<" | "<<xyz[2]<<" | \n";
			
			int row= EN_id(vertex);

			// search entity in which the point above are
			pEntity element = (pEntity)Octree_Search(xyz,theOctree);
			if (!element){
				cout << "Entity not found! Exiting...\n";
				exit(1);
			}

			//get coordinates
			vertice_1 = element->get(0,0);
			V_coord(vertice_1,ptn1);
			int id1 = EN_id(vertice_1);
			
			//cout<<"  Node "<<id1<<" : "<<ptn1[0]<<" | "<<ptn1[1]<<" | "<<ptn1[2]<<" | \n";
			
			vertice_2 = element->get(0,1); V_coord(vertice_2,ptn2); int id2 = EN_id(vertice_2);
			
			//cout<<"  Node "<<id2<<" : "<<ptn2[0]<<" | "<<ptn2[1]<<" | "<<ptn2[2]<<" | \n";
			
			vertice_3 = element->get(0,2); V_coord(vertice_3,ptn3); int id3 = EN_id(vertice_3);
			
			//cout<<"  Node "<<id3<<" : "<<ptn3[0]<<" | "<<ptn3[1]<<" | "<<ptn3[2]<<" | \n";

			//calculate areas
			double a1 = F_area(xyz,ptn2,ptn3);
			
			//cout<<"Area1: "<<a1<<endl;
			
			double a2 = F_area(xyz,ptn1,ptn3);
			
			//cout<<"Area2: "<<a2<<endl;
			
			double a3 = F_area(xyz,ptn1,ptn2);
			
			//cout<<"Area3: "<<a3<<endl;
			
			double a = a1 + a2 + a3;
			
			//cout<<"Area: "<<a<<endl;
			
			a1=a1/a; a2=a2/a; a3=a3/a;
			
			//cout<<"      Coeff1: "<<a1<<endl;
			//cout<<"      Coeff2: "<<a2<<endl;
			//cout<<"      Coeff3: "<<a3<<endl;

			//calculate coefficients
			setGeomCoeff(row,0,a1);
			setGeomCoeff(row,1,a2);
			setGeomCoeff(row,2,a3);
		}

		VIter_delete(vit);
	}
	else {

		//variables
		pEntity vertice_1,vertice_2,vertice_3,vertice_4; //element vertices
		double ptn1[3],ptn2[3],ptn3[3],ptn4[3];

		//Loop over vertices
		VIter vit = M_vertexIter(mesh);
		while ( pEntity vertex = VIter_next(vit) ){

			// get coordinate of a point from the mesh to be interpolated
			double xyz[3];
			V_coord(vertex,xyz);
			int row= EN_id(vertex);
			// search entity in which the point above are
			pEntity element = (pEntity)Octree_Search(xyz,theOctree);
			if (!element){
				cout << "Entity not found! Exiting...\n";
				exit(1);
			}

			//get coordinates
			vertice_1= element->get(0,0); V_coord(vertice_1,ptn1);
			vertice_2= element->get(0,1); V_coord(vertice_2,ptn2);
			vertice_3= element->get(0,2); V_coord(vertice_3,ptn3);
			vertice_4= element->get(0,3); V_coord(vertice_4,ptn4);

			//calculate areas
			double v1 = R_Volume(xyz,ptn2,ptn3,ptn4);
			double v2 = R_Volume(xyz,ptn1,ptn3,ptn4);
			double v3 = R_Volume(xyz,ptn1,ptn2,ptn4);
			double v4 = R_Volume(xyz,ptn1,ptn2,ptn3);
			double v = v1 + v2 + v3 + v4;
			v1= v1/v; v2= v2/v; v3= v3/v; v4= v4/v;
			//calculate coefficients
			setGeomCoeff(row,0,v1);
			setGeomCoeff(row,1,v2);
			setGeomCoeff(row,2,v3);
			setGeomCoeff(row,3,v4);
		}

		VIter_delete(vit);
	}
}


void DataTransfer::calculate_LinearInterpolation(pMesh mesh,pMesh mesh1, int dim){

	//variables
	int id, row;
	double NodeValue, GeomCoeff, linear;
	pEntity v;
	
	//Loop over vertices (from the mesh to where the data is being transfered)
	VIter vit = M_vertexIter(mesh);
	while ( pEntity vertex = VIter_next(vit) ){

		//get the node coordinate 
		double xyz[3];
		V_coord(vertex,xyz);
		//cout<<xyz[0]<<" | "<<xyz[1]<<" | "<<xyz[2]<<" | \n";

		//search entity in which the node is
		pEntity element = (pEntity)Octree_Search(xyz,theOctree);
		if (!element){
			cout << "WARNING:  Entity not found!\n";
			//exit(1);
		}

		//interpolate     
		//get row that contains the coefficients to be used in the interpolation
		row= EN_id(vertex);
		linear=0;
		for(int i=0;i<(dim+1);i++){ 
		
			double p[3];
			
			//get node id 			
			v = element->get(0,i);
			int id= EN_id(v);
			//V_coord(v,p);
			
			//cout<<" No "<<id<<" : "<<p[1]<<" | "<<p[2]<<" | "<<p[3]<<endl;
			
			NodeValue= getNodeValue(id);
			
			//cout<<"   Scalar: "<<NodeValue<<endl;
			
			GeomCoeff= getGeomCoeff(row,i);
			linear+=NodeValue*GeomCoeff;
	
			//cout<<"   LINEAR + "<<NodeValue<<" * "<<GeomCoeff<<" = "<<linear<<"\n"<<endl;    
		} 

		setInterpolatedValues(row,linear);
		//EN_attachDataDbl(vertex, MD_lookupMeshDataId("InterpLinear"),linear);
	}
	VIter_delete(vit);
}



double DataTransfer::calculate_QuadraticInterpolation(pMesh theMesh, int dim){
  
  //cout << __LINE__ <<":\t" << __FILE__ << endl;
  //variables
	int id, row;
	double NodeValue, GeomCoeff, residual, coord[3], linear, quadratic;
	pEntity v;

	//Loop over vertices (from the mesh to where the data is being transfered)
	VIter vit = M_vertexIter(theMesh);
	while ( pEntity vertex = VIter_next(vit) ){
			  
		//get the node coordinate 
		double xyz[3];
		V_coord(vertex,xyz);
		row = EN_id(vertex);
		
		//cout << __LINE__ <<":\t" << __FILE__ << endl;

		//search entity in which the node is
		pEntity element = (pEntity)Octree_Search(xyz,theOctree);
		if (!element){
			cout << "WARNING:  Entity not found!\n";
			//exit(1);
		}
		
		//calculate residual
		residual = 0;
		for( int i = 0; i < dim+1; i++){
		  
		 // cout << __LINE__ <<":\t" << __FILE__ << endl;
		  
		  v = element->get(0,i);
		  int id= EN_id(v);
		  V_coord(v,coord);
		  
		  //cout << __LINE__ <<":\t" << __FILE__ << endl;
		  
		  GeomCoeff= getGeomCoeff(row,i);
		  double dx = getGrad(id,1);
		  double dy = getGrad(id,2);
		  double dx2 = getGrad(id,3);
		  double dy2 = getGrad(id,4);
		  double dxdy = getGrad(id,5);
		
		  double Res = (xyz[0]-coord[0])*dx + (xyz[1]-coord[1])*dy + (xyz[0]-coord[0])*(xyz[0]-coord[0])*dx2/2 
		  
				+ (xyz[1]-coord[1])*(xyz[1]-coord[1])*dy2/2 + (xyz[0]-coord[0])*(xyz[1]-coord[1])*dxdy;
			    
		  residual += GeomCoeff*Res;  
		  }
		  
		  //cout<<"Residual: "<<residual<<endl;
		  //calculate quadratic interpolation
		  linear = getInterpolatedValues(row);
		  quadratic = linear + residual;
		  setInterpolatedValues(row,quadratic);
		  
	}
	VIter_delete(vit);
	
	//cout << __LINE__ <<":\t" << __FILE__ << endl;
}
 

void DataTransfer::showScalarField() const{
	//for (int i=0; i<)
}
